Benchmarking Iterative Projection Algorithms for Phase Retrieval 
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Iterative projection algorithms for phase retrieval are tested on two simple 'toy' models. The 
result provides useful insights in the behavior of these algorithms. 



Iterative transform methods pioneered by Gerchberg 
and Saxton y|, are well established techniques for iter- 
atively recovering the phase from the knowledge of the 
diffraction amplitude. The development of iterative al- 
gorithms with feedback in the early nineteen-eighties by 
Fienup produced a remarkably successful optimization 
method capable of extracting phase information from ad- 
equately sampled intensity data 00- [1 Finally, the im- 
portant theoretical insight that these iterations may be 
viewed as projections in Hilbert space 0, @ has allowed 
theoreticians to analyze and improve on the basic Fienup 
algorithm 0, IE 0, 

These algorithms try to find the intersection between 
two sets, typically the set of all the possible objects with 
a given diffraction pattern (modulus), and the set of all 
the objects that are constrained within a given area called 
support (or solvent in crystallography). The search for 
the intersection is based on the information obtained by 
'projecting' the current estimate on the two sets. An 
error metric exists to characterize the distance between 
the current estimate and a given feasibility set. The er- 
ror metric and its gradient are used in conjugate gradient 
(CG) based methods such as SPEDEN jlj. A projector 
P is an operator that takes to the closest point of a set 
from the current point p. A repetition of the same projec- 
tion is equal to one projection alone (P 2 = P). Another 
operator used here is the reflector R = 2P — I. We con- 
sider two sets, S (support) and M (modulus). The sup- 
port constraint is convex, while the modulus constraint 
is non-convex. Problems arise for non-convex sets, where 
projections become multivalued [l2j . 

The support projector P s acts on the object density p 
by setting to the density of the object outside a given 
region. The modulus projector P m acts on the density p 
in the Fourier domain p by forcing the modulus \p\ to be 
equal to the known one m, but keeping the phase of the 
current object in the Fourier domain p. This operator 
is demonstrated to be a projector on the non-covex set 
of the magnitude constraint [l2^. The same paper dis- 
cusses the problems of multi-valued projections for non- 
convex sets, which do not statisfy the requirements for 
gradient-based minimization algorithms, and the related 
nonsmoothness of the squared set distance metric, which 
may lead to numerical instabilities. See also [13] for a 
follow-up discussion on the non-smooth analysis. 

Several algorithms based on these concepts have now 
been proposed and a visual representation of their be- 
haviour is usefull to characterize the algorithm in various 
situations, in order to help chose the most appropriate 



one for a particular problem. 

The following algorithms require a starting point p°, 
which is generated by assigning a random phase to the 
measured object amplitude in the Fourier domain \p\. 
The first algorithm called Error Reduction (ER) (Ger- 
chberg and Saxton fi|) (see also Alternating Projections 
Onto Convex Sets [l4[ or Alternating Projections Onto 
Nononvex Sets p|) is simply: 
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by projecting back and forth between two sets, it con- 
verges to the local minimum (gradient type). The eigen- 
values of the support projectors are and 1, with cor- 
risponding eigenvectors the pixels outside and inside the 
support. Replacing the support projector P s with its 
reflector R s = 2P S — I, the corrisponding eigenvalues 
become -1 and 1, i.e. the charge density p outside the 
support is multiplied by -1. This algorithm is called sol- 
vent flipping in crystallography |l5j | : 
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It is often used in conjunction of the ER algorithm, al- 
ternating several HIO iterations and one ER iteration 
(HIO(20)+ER(1) in our case). Difference Map with 
7i = — /3 _1 , 72 = 0, which requires 4 projections 
(two time-consuming modulus constraint projections): 
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The Averaged Successive Reflections (ASR) [8j is: 
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The Hybrid Projection Reflection (HPR) 9] is derived 
from a relaxation of ASR: 
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It is equivalent to HIO if positivity is not enforced but it 
is written in a recursive form, instead of a case by case 
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FIG. 1: The basic features of the iterative projection algorithms can understood by this simple model of two lines intersecting 
jl(a)| . The aim is to find the intersection. The ER algorithm and the Solvent flipping algorithms converge in some gradient 
type fashion (the distance to the two sets never increases), with the solvent flip method being slightly faster when the angle 
between the two lines is small. HIO and variants move slightly in the direction where the gap between the two proje ctions 
decreases, but at the same time in the direction of the gap, following a spiral path. When the two lines do not intersect | |l(b)| 
HIO and variants keep moving in the direction of the gap. ER, Solvent Flipping and RAAR converge at (or close to) the local 
minimum. 



form such as Eq. [21 Finally Relaxed Averaged Alternating 
Reflectors RAAR (previously named RARS) [l(| 
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For j3 = 1, HIO, HPR, ASR and RAAS coincide. 

The first test is performed on the simplest possible 
case: find the intersection between two lines. Fig. 
shows the behavior of the various algorithms, The two 
sets are represented by a horizontal blue line (support) 
and a tilted black line (modulus). ER simply projects 
back and forth between these two lines, and moves along 
the support line in the direction of the intersection. Sol- 
vent Flip projects onto the modulus, 'reflects' on the sup- 
port, and moves along the reflection of the modulus con- 
straint onto the support. The solvent flipping algorithm 
is slightly faster than ER due to the increase in the an- 
gle of the projections and reflections. HIO and variants 
(ASR, Difference Map, HPR and RAAS) move in a spiral 
around the intersection eventually reaching the intersec- 
tion. For similar (3 RAAS behaves somewhere in between 
ER and HIO with a sharper spiral, reaching the solution 
much earlier. Alternating 20 iteration of HIO and 1 of ER 
(HIO(20)+ER(1)) considerably speeds up convergence. 

W hen a gap is introduced between the two lines (Fig. 
1(b) I so that the two lines do not intersect, HIO and vari- 



ants move away from this local minimum in search for 
another 'attractor' or local minimum. This shows how 



these algorithms escape from local minima and explore 
the HUbert space for other minima. ER, Solvent Flip, 
RAAS converges to or near the local minimum. By vary- 
ing (3 RAAS becomes a local minimizer for small /?, and 
becomes like HIO for /3 ~ 1. ER, solvent flip HIO+ER 
converge to the local minimum. The properties of SPE- 
DEN cannot be fully apreciated in these examples since 
the support constraint is represented by a one dimen- 
sional set, and this conjugate gradient method is designed 
for multidimensional minimization. However it is impor- 
tant to remark that such algorithm converges quadrati- 
cally to the local minimum, reaching the intersection in 
a single step when it exists, and the local minimum when 
a gap is introduced between the two sets in fewer steps 
than any of the other algorithm described here. 

A more realistic example is shown in Fig. [21 Here the 
circumference of two circles represent the modulus con- 
straint, while the support constraint is represented by a 
line. The two circles are used to represent a non-convex 
set with a local minimum. It is difficult to represent a 
true modulus constraint in real space. For a represen- 
tation of the modulus constraint in reciprocal space see 
|12j . The advantage of this example is the simplicity in 
the 'modulus' projector operator (it projects onto the 
closest circle). Although a real modulus constraint pro- 
jector is not as simple as the one used in this example, 
there are similarities: each Fourier space point provides 
an n-dimensional ellipsoid type equation. 
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We start from a position near the local minimum. ER, 
solvent flip and HIO+ER all fall into this trap (Fig. 
2(a) I, although increasing the interval between ER itera- 
tions in the HIO+ER algorithm would allow it to escape 
this local minimum. HIO and variants move away from 
the local minimum, 'find' the other circle, but converge 
to the center of the circle, with all but Diff. Map. not 
reaching a solution. In the center of the circle the pro- 
jection on the modulus constraint becomes 'multivalued', 
and its distance metric is 'nonsmooth'. The introduction 
of a small a random number added to the resulting solu- 
tion at every step allows all the HfO-typ e cod es to escape 
stagnation and find the solution (Fig. 2(b) I. The ran- 
dom number can be as low as the numerical precision 
of the computer. For (3 reduced to .9, RAAS would not 
reach the solution, but converge close to the local mini- 
mum. As a latest test in this series Fig. |2(c)| shows the 
behaviour of the algorithms when the support is tangent 
to the circle, the two solutions coincide, and the the two 
constraints are parallel. The only algorithm to reach the 
solution is RAAS, but HIO+ER would also reach the so- 
lution if the interval between ER steps was sufficiently 
large. 



I. POSITIVITY 

The situation changes slightly when we consider the 
positivity constraint. The previous definitions of the al- 
gorithms still apply just replacing P$ with Ps+- 
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otherwise. 
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The only difference is HIO which becomes: 
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if x £ S and P r 
otherwise. 



Fig. 3(a) shows that HIO bouches at the x 



(9) 
axis. 

As the positivity constraint gets closer to the solution, 
none of the algorithms converges to the solution (Fig. 
3(b) I, with the HIO- type algorithms bouncing between 
the regions closer to the two ci rcles. Only Difference 
Map for (3 > I converges (Fig. |3(c)|). Also HIO+ER 



would reach the solution for larger intervals between ER 
iterations. 



II. CONCLUSIONS 

ER is a simple but powerfull local minimizer, HIO and 
variants are very powerfull in escaping local minima, but 
in several situations fail to converge. When positivity 
is introduced, the recursive version of HIO (HPR) con- 
verges more 'smoothly' to the solution without bounch- 
ing on the x = axis. Alternating between HIO and 



ER with the correct intervals would have worked in all 
the examples shown above. RAAS is a good (single pa- 
rameter) way to change from 'global' to local minimizer, 
although it seems better to start from a high value of 
(3 and decrease it afterwards. Difference Map is succes- 
full in a few more of the examples shown above for the 
proper choiche of /3, however it involves 2 time consuming 
modulus constraint operations. To find when stagnation 
occours one can monitor the distances (using the proper 
error metrics) of the current solution before and after 
aplying various projectors, or monitor the autocorrela- 
tion between two succesive reconstructions. SPEDEN, a 
conjugate gradient based method, reaches a local mini- 
mum with quadradic convergence, and provides such kind 
of information. 

The Solvent flipping algorithm does not show much 
success in the examples shown above. Despite this it was 
used to improve images |l5j . and in a modified form to 
solve 3D structures ab-initio 01 • Perhaps the reason 
for the latter success has more to do with the thresh- 
old constraint used. A threshold projector multiplies by 
everything that is below a given threshold, while a 
threshold reflector multiplies it by -I. The application 
of this constraint has been proven succesfull in obtaining 
ab-initio solution in many circumstances applied in a va- 
riety of ways. Apart from Charge Flipping, which uses a 
threshold reflector fljl in electron density modification 
procedure with SIR |18(. in atomicity constraint |l6j| . a 
modified histogram constraint using the Difference Map 
(H. He, private comm.), and for complex valued objects, 
in the shrink-wrap algorithm using an updated support 
constraint by thresholding the current object reconstruc- 
tion at low resolution [l9| . 

By compressing the image in small spots and by 
increasing the flat region, these algorithms based on 
the threshold constraint resemble the maximum entropy 
method |20j , which tries to reduce the amout of informa- 
tion in an image and maximize the number of pixels of 
equal value (solvent). 

What distinguishes shrink-wrap and SIR algorithms 
is that they somehow solve the low resolution first, and 
gradually introduce high resolution information while 
slowly updating the low resolution information. Ficnup 
has also shown that starting from a low resolution im- 
age, slowly increasing resolution improves the algorithm 
|2ll | . Also SPEDEN starting from a low resolution target 
has been shown to slowly extend the correct phase in- 
formation to higher resolution. One possible explanation 
could be that the set based on low resolution is 'more 
smooth'. Perhaps also slowly increasing the gray lev- 
els in an image, or the number of possible phases of its 
Fourier transform could improve convergence by gradu- 
ally introducing more degrees of freedom in the resulting 
image, although a set defined by 'quantized' levels is also 
non-convex rendering it counterproductive. 
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FIG. 2: The horizontal line represents a support constraint, while the two circles represent a non convex constraint, i.e. the 
modulus constraint. The dashed line divides the region closer to one circle from the other. The starting point is on the circle 
to the right, possessing a local minimum distance to the line, (a) The gradient-type (ER and Solvent Flip) algorithms converge 
to the local minimum, while HIO and variants move away from the local minimum in the direction of the gap (vertical) untill 
they reach the region where the second circle is closer (delimited by the dashed line) . From here they try to move in the same 
spiral-like path of the two lines (Fig. untill they reach the point where the projecton on the circle and the line are parallel, 
and start moving toward the the center of the circle which has the correct solution. They stagnate in the center of the circle 
where the projection is multivalued. Only the Diff. Map reaches one of the two solutions. The addition of a small value of 
the order of the numerical precision after each iteration solves this stagnation (b). When one of the circles just touches the 
other constraint most algorithms either get stuck near the local minimum or stagnate. RAAS is the only one that reaches the 
vicinity of the solution (c). 
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(c) 

FIG. 3: (a) The starting point is again on the circle to the right, close to a local minimum. HIO and variant move away from 
the local minimum in the direction of the gap untill they reach the region where the circle to the left is closer. Instead of moving 
in a spiral like fashon, the iterations move close to the dotted line joining the center of the left circle to the origin, except for 
HIO that bounches on the x=0 axis, (b) The solution is very close to 0, and the dotted line originating from the circle to the 
left and passing by the origin becomes more tilted. The various algorithms after moving in the vertical direction away from the 
local minimum, reach the dashed line and start moving toward the tilted dotted line, falling back in the region closer to the 
first minimum. Tthese algorithms bounce between the regions closer to each circle without reaching the solution. With > 1, 
i.e. inverting the order of the operators, Diff. Map converges, and RAAS diverges, while HIO HPR and ASR stagnate. 
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